Suitability changes of Citrus medica L. var. sarcodactylis Swingle, a medicine-food plants affected by climate warming using the optimized MaxEnt model

Climatic variables are important conditions for plant growth, development and reproduction. Citrus medica L. var. sarcodactylis Swingle (Rutaceae: Citrus) is one of the traditional bulk Chinese medicinal materials in China with the effects of bacteriostasis, anti-inflammatory, anti-oxidation, anti-cancer cells, regulating the immun. Analyzing the impact of climate change on geographical distribution of C. medica L. var. sarcodactylis can provide strong support for its production layout and agricultural zoning. In our paper, MaxEnt and ArcGIS were applied to simulate the suitable areas of C. medica L. var. sarcodactylis in China from the perspectives of bioclimate, soil, topographic factors and human activities, and the future climate scenarios generated by global climate models (GCMs) were selected to predict its suitable areas in 2050s and 2090s. Results showed that, 1) Under current climate condition, areas of the total, most, moderately and poorly suitable habitats of C. medica L. var. sarcodactylis in China were 177.36×104 km2, 22.27×104 km2, 51.96×104 km2 and 103.13×104 km2 respectively. The range of the most suitable habitat was the narrowest, which was located in the middle east of Sichuan, western Chongqing in the upstream of the Yangtze River Basin, southern Guizhou and western Guangxi in the upstream of the Pearl River Basin, central and southern Yunnan and Southeast Tibet in the Middle-Lower reaches of the Southwest River Basin and western Taiwan. 2) Under the future climate change scenarios, the total suitable area showed a significant increase trend in 2090s, and the change of most, moderately and poorly suitable habitats showed no obvious law. 3) Under SSP1-2.6, SSP2-4.5 and SSP5-8.5 scenarios, the centroid of the most suitable habitat of C. medica L. var. sarcodactylis would move to the northwest, southeast and southwest respectively.


Introduction
Climatic variables are conditions for plant growth, development and reproduction [1,2]. The distribution of plants depends largely on abiotic factors, such as precipitation, temperature, soil and altitude [3,4]. Studying the impact of climate change on plant distribution is helpful to agricultural production and biodiversity protection, and can promote the sustainability of ecosystem. The drastic change of climate will lead to the change of the original habitat of plants, and then affect their suitable range [5][6][7]. Scientific evidences show that climate events such as climate warming, precipitation pattern change and atmospheric CO 2 concentration rise have made a significant influence on species diversity, genetic diversity and landscape diversity [8,9]. The formation of traditional Chinese medicine is closely related to climate, and climate change will affect its suitable growth areas [10,11]. Among the nearly 200 genuine medicinal materials in China, some genuine areas have remained stable, such as Chaenomeles speciosa and Fritillaria cirrhosa, while some have undergone corresponding changes, such as Alisma orientale and Citrus aurtantium [10].
The global climate change characterized by the increase of temperature and the change of precipitation makes the natural system, biological system and human health system oscillate [12][13][14]. Agricultural production, which relies heavily on natural resources, shows obvious vulnerability under climate change. This is because climate determines the production potential and yield of crops to a certain extent [2]. Temperature rise will increase the heat stress of most regions and crops in the world, and the instability of water stress due to precipitation fluctuation will pose a threat to crop production [15,16]. The rainstorm, flood, high temperature and drought disasters brought by climate change have a direct impact on agricultural production and pose a great challenge to people's growing demand and desire for a better life [17,18]. It is of great significance for agricultural production, social production and lifestyle to analyze the spatiotemporal change trend and distribution characteristics of temperature and precipitation under the background of climate change [19,20]. Teixeira et al. [21] found through simulation that under the influence of changes in heat stress under the A1B emission scenario, the area of the world's crop suitable areas will change from 2071 to 2100, and the change in the area of corn and wheat will be much larger than that of rice and soybean. Tang et al. [22] calculated the spatial-temporal distribution of reference crop evapotranspiration in the Huang-Huai-Hai Plain under the main climate scenarios in the future by using temperature and precipitation data. The results provide basic data support for scientific allocation of agricultural water resources and scientific response to the impact of climate change on agricultural production.
Citrus medica L. var. sarcodactylis Swingle (Rutaceae: Citrus) is one of the traditional bulk Chinese medicinal materials in China, which was first contained in the first TCM treatise Shengnong Bencao Jing (also Shennon, g's Classic of Materia Medica) [23]. The roots, stems, leaves, flowers and fruits of C. medica L. var. sarcodactylis can be used as medicine, and have unique functions of curing and preventing diseases, nourishing health and prolonging life [24]. Modern scientific experiments have shown that C. medica L. var. sarcodactylis has the effects of bacteriostasis, anti-inflammatory, anti-oxidation, anti-cancer cells, regulating the immune system [25]. At present, it is mainly used in medicine, food and daily chemical industries [26]. C. medica L. var. sarcodactylis is native to India and widely cultivated in France, Italy, Germany, the United States and Southeast Asia [23]. The planting areas in China are mainly distributed in Guangdong, Fujian, Chongqing, Sichuan and Zhejiang, especially Jiangjin in Chongqing and Zhaoqing in Guangdong, which have the largest planting area and the highest yield [27].
In recent years, the maximum entropy model (MaxEnt) has been applied to the study of habitat zoning of various medicinal plants under climate change scenarios since its good prediction ability. She et al. used MaxEnt to simulate the suitable distribution of the important medicinal resource plant Notopterygium incisum in the Three Rivers Headwater Region of China in the present and 2050s, and believed that the results would be conducive to the scientific protection and rational utilization of its resources [28]. Zhang et al. selected MaxEnt to compare the changes of the suitable habitat of wild Anemarrhena asphodeloides in China under the future climate change scenarios, and providing a theoretical basis for the protection of this medicinal plant and the scientific and reasonable increase of planting area [29]. The above research proves that MaxEnt model can be well applied to the study of habitat suitability analysis of medicinal plants.
With the increase of government investment in the traditional Chinese medicine industry, the cultivation area of C. medica L. var. sarcodactylis has been expanding [23]. However, without the guidance of scientific planting zoning and development planning, farmers planted it in large areas in non optimal growth areas, resulting in a series of problem. Therefore, it is necessary to analyze the suitability of C. medica L. var. sarcodactylis in China. In this paper, MaxEnt and ArcGIS were used to simulate the suitable areas of C. medica L. var. sarcodactylis in China from the perspectives of bioclimate, soil, topographic factors and human activities, and the climate change scenarios generated by global climate models (GCMs) were selected to predict its suitable areas in 2050s (2041-2060) and 2090s (2081-2100). The main purposes of this study are: 1) to provide scientific guidance for the reasonable planning and high-quality cultivation of C. medica L. var. sarcodactylis in China. 3) to master the impact of climate change on the planting suitability of C. medica L. var. sarcodactylis.

Occurrence data of C. medica L. var. sarcodactylis
Occurrence data of C. medica L. var. sarcodactylis were acquired from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), the Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn/), field survey, and literature. Referring to the methods in literature [30,31], we processed the distribution records of C. medica L. var. sarcodactylis. Firstly, Baidu coordinate picking system (https://api.map.baidu.com/lbsapi/getpoint) was used to determine the longitude and latitude of records accurate to the town level. Second, ENMTOOL was applied to process the data to ensure the uniqueness of distribution points in each grid. After the above procedures, 123 records of C. medica L. var. sarcodactylis were retained for the establishment of MaxEnt (Fig 1).

Environmental data
In this study, 27 initial environment variables with a resolution of 2.5 arcminutes were selected (S1 Table). BCC-CSM2-MR (Beijing Climate Center Climate System Model) developed by China National Climate Center in CMIP6 which has been proved to be more suitable for China's climate change characteristics was selected as the future climate model [32]. SSP1-2.6 (the lowly greenhouse gas emission scenario), SSP2-4.5 (the moderately greenhouse gas emission scenario) and SSP5-8.5 (the highly greenhouse gas emission scenario) in the Shared Socio-economic Pathways (SSPs) scenarios that can more scientifically describe future climate change were selected [33]. In order to reduce the influence of multi-collinearity among 19 bioclimatic variables, Pearson correlation analysis method was used [31,34]. Firstly, MaxEnt was used to calculate the percent contribution of 19 environmental variables, and the variables whose percent contribution rate was greater than 0 were retained (S1 Table). Thereafter, the Pearson's coefficients between two variables with percent contribution greater than 0 corresponding to 123 occurrence data of C. medica L. var. sarcodactylis were analyzed using SPSS (S2 Table).
Thirdly, by comparing the percentage contribution of the variables with the absolute value of the coefficient greater than 0.85, the higher one was retained. Finally, in addition to elevation, a total of 17 variables were selected to establish the prediction model of C. medica L. var. sarcodactylis (S1 Table).

Modelling process
MaxEnt software (Version 3.4.4, developed by Phillips et al. [39]) operation procedure was as follows: 1) The occurrence data of C. medica L. var. sarcodactylis in "CSV" format and the environmental variable in "ASC" format were imported into the "sample" and "environmental layers" data boxes respectively. 2) "Create response curves" and "Do jackknife to measure variable importance" were selected respectively to analyze the relationship between variables and presence probability of C. medica L. var. sarcodactylis and measure the importance of variables. 3) In the initial model, "Random test percentage" was set to 25%, while in the reconstructed model, "random seed" was selected, and the "replicates" was set to 10 [31,34,35]. 4) The Kuenm software package of R language was used to optimize the regularization multiplier (40, with values of 0.1~40) and feature combination (L, Q, P, T, H, LQ, LP, LT, LH, QP, QT, QH, PT, PH, TH, LQP, LQT, LQH, LPT, LPH, QPT, QPT, QPH, QTH, QTH, LQPT, LQPH, LQPH, LQTH, LQTH and LQPTH) of the model, and the optimal setting of the minimum information criterion AICc value (delta.AICc) among 1160 results was selected [36]. When the continuous prediction results are converted into the Boolean values of the "suitable habitat" and "unsuitable habitat", it is critical to select the appropriate threshold. In this study, the probability value P corresponding to the maximum sum of sensitivity and specificity was taken as the threshold value, that was, the probability of species existence p � P was taken as the suitable hbitat, and p<P was taken as the unsuitable habitat. The suitable habitat was subdivided according to the following standards: poorl suitable habitat (P�p<0.33), moderately suitable habitat (0.33�p< 0.66) and most suitable habitat (p�0.66) [6,37,38].

Suitable habitat under current climate condition
After reclassification and format conversion of the results output by MaxEnt, the area calculated by ArcGIS showed that the total suitable area of C. medica L. var. sarcodactylis in China was 177.36×10 4 km 2 , accounting for 18.47% of the total area of the country, of which the distribution range of the most suitable habitat was the narrowest, with an area of 22.27×10 4 km 2 , accounting for 12.55% of the total suitable habitat. River Basin and central Guangdong in the lower reaches, southern Yunnan in the lower reaches of the Southwest River basin and most of Hainan (Fig 2). The poorly suitable habitat of C. medica L. var. sarcodactylis were mainly distributed in the northeast and south of Sichuan, the north of Chongqing, the south of Shaanxi, the west of Hubei, the west of Hunan, the northeast and northwest of Guizhou in the middle and upper reaches of the Yangtze River Basin, the east of Jiangxi in the lower reaches of the Yangtze River, most of Fujian in the southeast River Basin, most of Guangxi in the upper reaches of the Pearl River Basin and Guangdong in the lower reaches, southeast of Tibet and north central Yunnan in the lower reaches of the Southwest River basin, covering an area of 103.13×10 4 km 2 , accounting for 58.15% of the total suitable habitat (Fig 2).
Under SSP5-8.5, the area of total suitable habitat was 199.08×10 4 km 2 (2050s) and 198.72×10 4 km 2 (2090s), which increased by 12.25% (2050s) and 12.04% (2090s) compared with current, respectively. The areas of most suitable habitats were 25.11×10 4 km 2 (2050s) and 23.56×10 4 km 2 (2090s), which increased by 12.79% (2050s) and 15.83% (2090s) compared with current, respectively. The areas of moderately suitable habitats were 56.84×10 4 km 2 (2050s) and 51.86×10 4 km 2 (2090s), which increased by 9.38% (2050s) and decreased by 0.2% (2090s) compared with current, respectively. The areas of poorly suitable habitats were 117.13×10 4 km 2 (2050s) and 123.3×10 4 km 2 (2090s), which increased by 13.57% (2050s) and 19.55% (2090s) compared with current, respectively (Fig 3E, 3F and 3G).  the suitability change of the suitable habitat in 2050s was obtained (Fig 4G). The results showed that eastern Sichuan, Chongqing, most of Guizhou, southern Yunnan, most of Guangxi, western Guangdong and Hainan were identified as stable habitats for the growth and distribution of C. medica L. var. Sarcodactylis (Fig 4G). The expand habitat was dotted on the map with an area of 2.24×10 4 km 2 , which was relatively concentrated only at the junction of Hubei, Chongqing, Hunan and Guizhou (Fig 4G). The area of the shrink habitat was 4.68×10 4 km 2 , which was more obvious on the map as the northeast of Fujian, the middle of Guangdong, and the junction of Guizhou, Guangxi and Yunnan (Fig 4G). Fig 4D-4F showed the stable, expand and shrink area of the suitable habitat in 2090s compared with current. The area proportion of the stable habitat to the current was 84.19% (SSP1-2.6), 85.37% (SSP2-4.5) and 82% (SSP5-8.5) respectively, the area proportion of the shrink habitat to the current was 15.81% (SSP1-2.6), 14.62% (SSP2-4.5) and 18% (SSP5-8.5) respectively, and the area proportion of the expand habitat to the future was 20.05% (SSP1-2.6), 18.74% (SSP2-4.5) and 19% (SSP5-8.5) respectively. By overlie Fig 4D-4F in ArcGIS, the suitability change of the suitable habitat in 2090s was obtained (Fig 5H). The results showed that the distribution pattern of the stable habitat was basically consistent with that of 2050s ( Fig 5H). Compared with 2050s, the expand habitat would increased significantly, especially at the junction of Sichuan, Yunnan and Guizhou. The shrink habitat was mainly located in northern Guangdong, southern Sichuan and northeastern Yunnan.   3˚N (2090s). From current to 2090s, the centroid would generally moved 57.89 km to the southeast (Fig 5). 3) Under SSP5-8.5 scenario, the centroid would move from 105.9˚E/26.68˚N (current) to the northeast by 19.23 km to 106.09˚E/26.7˚N (2050s), and then to the southwest by 35 .1 km to 105.74˚E/26.66˚N (2090s). From current to 2090s, the centroid would generally moved 15.88 km to the southwest (Fig 5). Table 1 showed that the AUC values of the training data and test data of the MaxEnt model under current situation were 0.969±0.001 and 0.952±0.013 respectively. The AUC values of the training data of the future climate scenario models were 0.963±0.002-0.969±0.002 and that of the test data were 0.945±0.014-0.955±0.018.

Discussion
Our simulation showed that the most suitable habitat of C. medica L. var. Sarcodactylis were distributed in the middle east of Sichuan, western Chongqing, southern Guizhou and western Guangxi, central and southern Yunnan and southeast Tibet, with an area of 22.27×10 4 km 2 . C. medica L. var. sarcodactylis is native to India, mainly cultivated in China, with few wild species. With the cultivation and introduction, C. medica L. var. sarcodactylis has been widely planted in Sihui, Chaoshan, Yunfu and Yunan of Guangdong, Tianlin, Longlin, Lingle, Guanyang, Daxin and Yongfu of Guangxi, Kunming, Yuxi, Chuxiong, Xinping, Yimen, Eshan and Pu'er of Yunnan, Jiangjin, Yongchuan, Yunyang and Kai counties of Chongqing, Anxian, Hejiang, Yibin, Muchuan, Ya'an, Hongya, Jiajiang, Qianwei and Xingjing of Sichuan [27]. The above areas have a high degree of overlap with our simulation, indicating the accuracy of MaxEnt. According to the data, the actual planting area of C. medica L. var. sarcodactylis is far smaller than its suitable area. MaxEnt only simulates the probability of existence of C. medica L. var. sarcodactylis under the combination of dominant natural environmental variables, and does not involve production management and the impact on fruit yield and quality. In the actual production, not only the dominant natural resources, but also the social and economic factors such as labor, field management, production cost and market demand should be considered [23,42], so that the actual production area is smaller than the suitable habitat.
Climate change is one of the most important threats to global biodiversity in this century. Climate change may lead to changes in the future geographical distribution pattern of species, exacerbate the reduction of biodiversity and the loss of germplasm resources, and even accelerate the extinction rate of species [43,44]. Therefore, it is of great significance for species protection and sustainable utilization of resources to carry out research on the distribution pattern and change of species suitable areas in the context of climate change. Studying the response of plant distribution pattern to climate change and understanding the relationship between plant climate demand and geographical distribution is of great significance to reveal the history of species formation, migration and diffusion, and to put forward reasonable utilization strategies and planting zoning [45][46][47]. Global climate change is an important driving factor to change the geographical distribution pattern of species [48]. Some species benefit from climate change, and their distribution will be expanded, while the habitats of some species on the contrary will be reduced [37,38,49,50]. Our research showed that under the SSP1-2.6 and SSP2-4.5 in 2050s, the total suitable area of C. medica L. var. sarcodactylis changed slightly, while under other scenarios, the area increased significantly. This indicated that the future climate change would be favorable to C. medica L. var. sarcodactylis as a whole, and this expanding trend had also been found in other medicinal plants. Xu et al. predicted the distribution of Thesium chinense in China, and their results revealed that the area of the highly suitable area would increased gradually and expanded from the Yangtze River basin to the Yellow River Basin [51]. Zhao et al. analyzed the geographical distribution of Tibetan medicine Lamiophlomis rotata in the Qinghai Tibet Plateau, and the results showed that the expanded area of the suitable habitat caused by climate change was about twice the degraded area, and the area in all provinces showed an increasing trend [13]. The higher precipitation under the high concentration emission scenario can reduce the restriction of precipitation factors on the distribution of species and expand the suitable habitat of species. On the contrary, the increased precipitation under the low concentration emission scenario can not reduce the limit, but will reduce the available water for species to absorb with global warming [12,14,52]. This is consistent with the change law of the total suitable habitat of C. medica L. var. sarcodactylis under climate change. We speculated that the reduction of the total suitable area of C. medica L. var. sarcodactylis in 2050s may be due to the excessive drought caused by the temperature rise. While with the increase of years, rainfall gradually accumulated, making some areas suitable. Trajectory analysis showed the centroid of the most suitable habitat of C. medica L. var. sarcodactylis would move to the northwest (SSP1-2.6), southeast (SSP2-4.5) and southwest (SSP5-8.5) respectively. This indicated that the impact of climate change on species distribution was uncertain, and the variation of temperature and precipitation under different scenarios may be the cause of this phenomenon. Xiang et al. found through simulation that under SSP1-2.6 scenario, the increasing trend of temperature rise rate and extreme precipitation rate in China from 2021 to 2100 was relatively flat, while under SSP2-4.5, SSP3-7.0 and SSP5-8.5 scenarios, they gradually increase with the increase of mode and time [53].
Since the 20th century, with the rapid growth of population and economy, human beings have had a great impact on their own living environment [54,55]. Climate change and human activities are considered to be the leading factors affecting the structure and function of terrestrial ecosystems [56,57]. Our analysis showed that human footprint was the key variables affecting the distribution of C. medica L. var. sarcodactylis, and the appropriate range was � 6.99. C. medica L. var. sarcodactylis is mostly cultivated in China, and the impact of agricultural activities on its distribution cannot be ignored. Cutting, grafting and high-pressure propagation are the main propagation methods of C. medica L. var. sarcodactylis [25]. In recent years, with the continuous optimization of agricultural industrialization, the continuous adjustment of agricultural planting structure, and the influence of market regulation, the planting area of C. medica L. var. sarcodactylis has expanded [58,59]. However, Yue et al. [60] investigated the resources of C. medica L. var. sarcodactylis in Guangdong, and the results showed that due to the influence of production technology, market conditions, economic benefits and other factors, the planting area of C. medica L. var. sarcodactylis in Guangdong was shrinking year by year, and even sporadic cultivation was found in only 13 of the 40 sample plots visited. The above results all prove the importance of human activities to the cultivation and distribution of C. medica L. var. sarcodactylis.

Conclusions
Under current climate condition, areas of the total, most, moderately and poorly suitable habitats of C. medica L. var. sarcodactylis in China were 177.36×10 4 km 2 , 22.27×10 4 km 2 , 51.96×10 4 km 2 and 103.13×10 4 km 2 respectively. The range of the most suitable habitat was the narrowest, which was located in the middle east of Sichuan, western Chongqing in the upstream of the Yangtze River Basin, southern Guizhou and western Guangxi in the upstream of the Pearl River Basin, central and southern Yunnan and Southeast Tibet in the Middle-Lower reaches of the Southwest River Basin and western Taiwan. Under the future climate change scenarios, the total suitable area showed a significant increase trend in 2090s, and the changes of most, moderately and poorly suitable habitats showed no obvious law. Under SSP1-2.6, SSP2-4.5 and SSP5-8.5 scenarios, the centroid of the most suitable habitat would move to the northwest, southeast and southwest respectively. This study predicted the response of suitable habitats of C. medica L. var. sarcodactylis to climate change, which is of great significance for its planting division and resource protection.
Supporting information S1